*Load up data file
use "`base'\Data Files and Code That Produced Them\P03 ED Analysis File", replace

*These coefficients come from "P11 Estimate population for weights.sas" and are the slope of a regression of population on age in years. So convert age in months in this file
*into years then multiply by slope to get weights. Add one so it is effectively proportion and deals with 0 issue. Stata will do the rescaling.  
gen pop_all_np =  1 + 0.00104*(months_21/12)
gen pop_m      =  1 + 0.00016718*(months_21/12)
gen pop_f_np   =  1 + 0.00196*(months_21/12)
 
 
local sufs = "all all_np f f_np m"  
local vars  = "visit  illness injury_or_alc  alcohol_any inj_by_self inj_by_oth inj_accident   "

*Create all fitted lines on either side of RD
foreach suf of local sufs {
   foreach var of local vars {  
     *Fit the left side polynomial;
      reg `var'_`suf'_r age_c age_c_sq [aweight=pop_`suf'] if age_c >= -2 & age_c < 0, robust 
      predict `var'_`suf'_left        if age_c >= -2 & age_c <= 0
     *Fit the right side polynomial;
      reg `var'_`suf'_r age_c age_c_sq [aweight=pop_`suf'] if age_c > 0 & age_c < 2, robust
      predict `var'_`suf'_right       if age_c >= 0 & age_c < 2
   }
}
  
  
keep if age_months >= 21 - 2  & age_months < 21 + 2   

sum alcohol_any_all_np_r inj_by_self_all_np_r inj_by_oth_all_np_r inj_accident_all_np_r

*Total split into two groups
gen dif = injury_or_alc_all_np_r - (alcohol_any_all_np_r + inj_by_self_all_np_r + inj_by_oth_all_np_r + inj_accident_all_np_r)
sum  dif
drop dif
tab inj_accident_all_np_r

#delimit ;
graph twoway  
 	        (scatter inj_accident_all_np_r age_months, mcolor(blue)  msymbol(O)  msize(.6) yaxis(1) 
			   yscale(titlegap(2)) 
			   ylabel(0(400)1200, axis(1) nogrid)) 
			(line inj_accident_all_np_left age_months, lwidth(thin) lcolor(blue) yaxis(1)) 
            (line inj_accident_all_np_right age_months, lwidth(thin) lcolor(blue) yaxis(1)) 
			  
	        (scatter alcohol_any_all_np_r age_months, mcolor(black)  msymbol(D)  msize(.4) yaxis(2) 
			   ylabel(0(50)150, axis(2) nogrid))
            (line alcohol_any_all_np_left age_months, lwidth(thin) lcolor(black) yaxis(2)) 
            (line alcohol_any_all_np_right age_months, lwidth(thin) lcolor(black) yaxis(2)) 

	        (scatter inj_by_self_all_np_r age_months, mcolor(red)  msymbol(Th)  msize(.8) yaxis(2)) 
            (line inj_by_self_all_np_left age_months, lwidth(thin) lcolor(red) yaxis(2)) 
            (line inj_by_self_all_np_right age_months, lwidth(thin) lcolor(red) yaxis(2)) 
			
	        (scatter inj_by_oth_all_np_r age_months, mcolor(green)  msymbol(Sh)  msize(.8) yaxis(2)) 
            (line inj_by_oth_all_np_left age_months, lwidth(thin) lcolor(green) yaxis(2)) 
            (line inj_by_oth_all_np_right age_months, lwidth(thin) lcolor(green) yaxis(2)) 			
             , 
			 title("Figure 2: Emergency Department Visits by Cause",size(medlarge)) 
			 xtitle("Age at Time of ED Visit",size(medlarge))  
			 ytitle("Accidental Injuries",axis(1) size(medlarge)) 
			 ytitle("Alcohol Intoxication & Deliberate Injuries", axis(2) size(medlarge)) 
			 legend(off) 
			 text(1150 20 "Accidental Injury",color(blue) size(medlarge)) 
			 text( 600 20 "Alcohol Intoxication", color(black) size(medlarge))
			 text( 290 20 "Deliberately Self Inflicted Injury", color(red) size(medlarge)) 
			 text( 920 20 "Deliberately Injured by Other Person", color(green) size(medlarge)) 
             graphregion(style(none) color(gs16))  
 
      ;
#delimit cr
grexportpdf "`base'\Code for Figures and Tables in Paper\Figure 2.pdf"



